1 Session Info

sessioninfo::package_info()
##  package        * version    date       lib source                         
##  ape              5.3        2019-03-17 [1] CRAN (R 3.6.0)                 
##  assertthat       0.2.1      2019-03-21 [1] CRAN (R 3.6.0)                 
##  backports        1.1.9      2020-08-24 [1] CRAN (R 3.6.2)                 
##  cli              2.0.2      2020-02-28 [1] CRAN (R 3.6.0)                 
##  coda             0.19-3     2019-07-05 [1] CRAN (R 3.6.0)                 
##  codetools        0.2-16     2018-12-24 [1] CRAN (R 3.6.1)                 
##  colorspace       1.4-1      2019-03-18 [1] CRAN (R 3.6.0)                 
##  crayon           1.3.4      2017-09-16 [1] CRAN (R 3.6.0)                 
##  DEoptimR         1.0-8      2016-11-19 [1] CRAN (R 3.6.0)                 
##  deSolve        * 1.27.1     2020-01-02 [1] CRAN (R 3.6.0)                 
##  digest           0.6.25     2020-02-23 [1] CRAN (R 3.6.0)                 
##  doParallel       1.0.15     2019-08-02 [1] CRAN (R 3.6.0)                 
##  dplyr            1.0.2      2020-08-18 [1] CRAN (R 3.6.2)                 
##  egonet         * 0.0.0.9000 2020-02-11 [1] Github (statnet/egonet@70d4aa2)
##  ellipsis         0.3.1      2020-05-15 [1] CRAN (R 3.6.2)                 
##  EpiModel       * 2.0.3      2020-11-09 [1] CRAN (R 3.6.2)                 
##  ergm           * 3.11.0     2020-10-14 [1] CRAN (R 3.6.2)                 
##  ergm.ego       * 0.6.1      2020-11-19 [1] CRAN (R 3.6.2)                 
##  evaluate         0.14       2019-05-28 [1] CRAN (R 3.6.0)                 
##  fansi            0.4.1      2020-01-08 [1] CRAN (R 3.6.0)                 
##  foreach          1.4.7      2019-07-27 [1] CRAN (R 3.6.0)                 
##  generics         0.0.2      2018-11-29 [1] CRAN (R 3.6.0)                 
##  ggplot2          3.3.2      2020-06-19 [1] CRAN (R 3.6.2)                 
##  glue             1.4.1      2020-05-13 [1] CRAN (R 3.6.2)                 
##  gtable           0.3.0      2019-03-25 [1] CRAN (R 3.6.0)                 
##  here           * 0.1        2017-05-28 [1] CRAN (R 3.6.0)                 
##  htmltools        0.5.0      2020-06-16 [1] CRAN (R 3.6.2)                 
##  iterators        1.0.12     2019-07-26 [1] CRAN (R 3.6.0)                 
##  knitr            1.29       2020-06-23 [1] CRAN (R 3.6.2)                 
##  lattice          0.20-38    2018-11-04 [1] CRAN (R 3.6.1)                 
##  lazyeval         0.2.2      2019-03-15 [1] CRAN (R 3.6.0)                 
##  lifecycle        0.2.0      2020-03-06 [1] CRAN (R 3.6.0)                 
##  lpSolve          5.6.13.3   2019-08-19 [1] CRAN (R 3.6.0)                 
##  magrittr         1.5        2014-11-22 [1] CRAN (R 3.6.0)                 
##  MASS             7.3-51.4   2019-03-31 [1] CRAN (R 3.6.1)                 
##  Matrix           1.2-17     2019-03-22 [1] CRAN (R 3.6.1)                 
##  munsell          0.5.0      2018-06-12 [1] CRAN (R 3.6.0)                 
##  network        * 1.16.1     2020-10-07 [1] CRAN (R 3.6.2)                 
##  networkDynamic * 0.10.0     2019-04-05 [1] CRAN (R 3.6.0)                 
##  nlme             3.1-140    2019-05-12 [1] CRAN (R 3.6.1)                 
##  pillar           1.4.6      2020-07-10 [1] CRAN (R 3.6.2)                 
##  pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 3.6.0)                 
##  purrr            0.3.4      2020-04-17 [1] CRAN (R 3.6.2)                 
##  R6               2.4.1      2019-11-12 [1] CRAN (R 3.6.0)                 
##  RColorBrewer     1.1-2      2014-12-07 [1] CRAN (R 3.6.0)                 
##  Rcpp             1.0.5      2020-07-06 [1] CRAN (R 3.6.2)                 
##  rlang            0.4.7      2020-07-09 [1] CRAN (R 3.6.2)                 
##  rle              0.9.2      2020-09-25 [1] CRAN (R 3.6.2)                 
##  rmarkdown        2.3        2020-06-18 [1] CRAN (R 3.6.2)                 
##  robustbase       0.93-5     2019-05-12 [1] CRAN (R 3.6.0)                 
##  rprojroot        1.3-2      2018-01-03 [1] CRAN (R 3.6.0)                 
##  scales           1.1.1      2020-05-11 [1] CRAN (R 3.6.2)                 
##  sessioninfo      1.1.1      2018-11-05 [1] CRAN (R 3.6.0)                 
##  statnet.common   4.4.1      2020-10-03 [1] CRAN (R 3.6.2)                 
##  stringi          1.4.6      2020-02-17 [1] CRAN (R 3.6.0)                 
##  stringr          1.4.0      2019-02-10 [1] CRAN (R 3.6.0)                 
##  tergm          * 3.6.1      2019-06-12 [1] CRAN (R 3.6.0)                 
##  tergmLite      * 2.2.1      2020-07-22 [1] CRAN (R 3.6.2)                 
##  tibble           3.0.3      2020-07-10 [1] CRAN (R 3.6.2)                 
##  tidyselect       1.1.0      2020-05-11 [1] CRAN (R 3.6.2)                 
##  trust            0.1-8      2020-01-10 [1] CRAN (R 3.6.0)                 
##  vctrs            0.3.2      2020-07-15 [1] CRAN (R 3.6.2)                 
##  withr            2.2.0      2020-04-20 [1] CRAN (R 3.6.2)                 
##  xfun             0.16       2020-07-24 [1] CRAN (R 3.6.2)                 
##  yaml             2.2.1      2020-02-01 [1] CRAN (R 3.6.0)                 
## 
## [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library

2 General Controls

need to update duration targets w/ estimated means from parametric model
but then will also need to use CKs correction in-simiulation find undelrying relationship length in presense of mortality

np <- detectCores()
# Set parallel.type
if (np==0) ptype=NULL else ptype='PSOCK'

dat <- setup_params()
ppopsize <- dat$ppopsize
MCMCinterval <- dat$mcmcInterval
MCMCburnin <- dat$mcmcBurnin

time.step <- dat$time.step
age.width <- dat$age.width
dissolution <- dat$dissolution
#mRate <- dat$mRate

duration.mar <- dat$duration.marriage
duration.cohab <- dat$duration.cohab
duration.other <- dat$duration.other

DissolutionMar <- dissolution_coefs(dissolution=dissolution, 
                                       duration=duration.mar)

DissolutionCohab <- dissolution_coefs(dissolution=dissolution, 
                                       duration=duration.cohab)

DissolutionOther <- dissolution_coefs(dissolution=dissolution, 
                                      duration=duration.other)

3 Some Adjustments

  • make age/sex attribute (i.e. 15F, 23M)
  • sample down to 10,000
ssize <- 10000
# sample down to 10,000
egodat_marriage <- sample(egodat_marriage, ssize)

# pull sampled egos, attribute updates, and update other alter lists (for consistent nodeset)
egos <- egodat_marriage$egos

egodat_cohab$alters <- egodat_cohab$alters[egodat_cohab$alters$ego %in% egos$ego,]
egodat_other$alters <- egodat_other$alters[egodat_other$alters$ego %in% egos$ego,]
egodat_once$alters <- egodat_once$alters[egodat_once$alters$ego %in% egos$ego,]

egos$ageF <- ifelse(egos$male==0, egos$age, 0)
egos$ageM <- ifelse(egos$male==1, egos$age, 0)

egodat_marriage$alters$ageF <- ifelse(egodat_marriage$alters$male==0, egodat_marriage$alters$age, 0)
egodat_marriage$alters$ageM <- ifelse(egodat_marriage$alters$male==1, egodat_marriage$alters$age, 0)

egodat_cohab$alters$ageF <- ifelse(egodat_cohab$alters$male==0, egodat_cohab$alters$age, 0)
egodat_cohab$alters$ageM <- ifelse(egodat_cohab$alters$male==1, egodat_cohab$alters$age, 0)

egodat_other$alters$ageF <- ifelse(egodat_other$alters$male==0, egodat_other$alters$age, 0)
egodat_other$alters$ageM <- ifelse(egodat_other$alters$male==1, egodat_other$alters$age, 0)

egodat_once$alters$ageF <- ifelse(egodat_once$alters$male==0, egodat_once$alters$age, 0)
egodat_once$alters$ageM <- ifelse(egodat_once$alters$male==1, egodat_once$alters$age, 0)

egodat_marriage$egos <- egos
egodat_cohab$egos <- egos
egodat_other$egos <- egos
egodat_once$egos <- egos


# grab binary deg dist for plots below, x2 for scaline to 20000 nodes
m <- with(egodat_marriage$egos[egodat_marriage$egos$male==1,], table(age, deg.mar))
m <- m[,2]*(ppopsize/ssize)
f <- with(egodat_marriage$egos[egodat_marriage$egos$male==0,], table(age, deg.mar))
f <- f[,2]*(ppopsize/ssize)

mc <- with(egodat_cohab$egos[egodat_cohab$egos$male==1,], table(age, deg.cohab))
mc <- mc[,2]*(ppopsize/ssize)
fc <- with(egodat_cohab$egos[egodat_cohab$egos$male==0,], table(age, deg.cohab))
fc <- fc[,2]*(ppopsize/ssize)

mo <- with(egodat_other$egos[egodat_other$egos$male==1,], table(age, deg.other))
mo <- mo[,2]*(ppopsize/ssize)
fo <- with(egodat_other$egos[egodat_other$egos$male==0,], table(age, deg.other))
fo <- fo[,2]*(ppopsize/ssize)

4 Marriage Network

4.1 Ergm.Ego Estimation & Setup

fit.marriage <- ergm.ego(egodat_marriage ~ edges + 
                         nodecov(~ageF) + 
                         nodecov(~ageF^2) +
                         nodecov(~ageF<30) + 
                         nodecov(~ageF>30) +
                         nodecov(~ageM) + 
                         nodecov(~ageM^2) +
                         nodecov(~ageM<30) + 
                         nodecov(~ageM>30) + 
                         absdiff(~sqrtage + dat$shift.mar*(male==0)) + 
                         nodefactor("deg.other.binary") +
                         offset(nodematch("male", diff = FALSE)) +
                         offset(nodefactor("olderpartnerM")) +
                         #offset(nodefactor("debuted", levels=1)) + 
                         offset("concurrent"), 
                         offset.coef = c(-Inf, -Inf, -Inf), #offset.coef = c(-Inf, -Inf, -Inf, -Inf),
                         control = 
                         control.ergm.ego(ppop.wt = "sample", ppopsize = ppopsize, 
                                          ergm.control = 
                                            control.ergm(MCMLE.maxit = 100,
                                                         MCMC.interval = MCMCinterval, 
                                                         MCMC.burnin = MCMCburnin,
                                                         parallel = np,
                                                         parallel.type = ptype)))

#nodefactor(~ageF, levels=5:9) + 
#nodefactor(~ageM, levels=7:10) + 

saveRDS(fit.marriage, here("fits", "fit.marriage.rds"))
fit.marriage <- readRDS(here("fits", "fit.marriage.rds"))
mpop <- environment(fit.marriage$ergm.formula)$popnw
marriage.netest <- ego.netest(fit.marriage, DissolutionMar)
summary(fit.marriage)
## Call:
## ergm(formula = ergm.formula, offset.coef = ergm.offset.coef, 
##     target.stats = m, eval.loglik = FALSE, control = control$ergm.control)
## 
## Iterations:  3 out of 100 
## 
## Monte Carlo MLE Results:
##                                           Estimate Std. Error MCMC % z value
## offset(netsize.adj)                      -9.903488   0.000000      0    -Inf
## edges                                   -18.415798   2.313703      0  -7.959
## nodecov.ageF                              0.658567   0.103593      0   6.357
## nodecov.ageF^2                           -0.008985   0.001596      0  -5.630
## nodecov.ageF<30                          -0.302788   0.503616      0  -0.601
## nodecov.ageF>30                          -0.696689   0.513842      0  -1.356
## nodecov.ageM                              0.581784   0.113492      0   5.126
## nodecov.ageM^2                           -0.007603   0.001686      0  -4.509
## nodecov.ageM<30                           0.381250   0.295776      0   1.289
## nodecov.ageM>30                           0.260988   0.294664      0   0.886
## absdiff.sqrtage+dat$shift.mar*(male==0)  -3.700792   0.141901      0 -26.080
## nodefactor.deg.other.binary.1            -4.394601   0.326862      0 -13.445
## offset(nodematch.male)                        -Inf   0.000000      0    -Inf
## offset(nodefactor.olderpartnerM.1)            -Inf   0.000000      0    -Inf
## offset(concurrent)                            -Inf   0.000000      0    -Inf
##                                         Pr(>|z|)    
## offset(netsize.adj)                       <1e-04 ***
## edges                                     <1e-04 ***
## nodecov.ageF                              <1e-04 ***
## nodecov.ageF^2                            <1e-04 ***
## nodecov.ageF<30                            0.548    
## nodecov.ageF>30                            0.175    
## nodecov.ageM                              <1e-04 ***
## nodecov.ageM^2                            <1e-04 ***
## nodecov.ageM<30                            0.197    
## nodecov.ageM>30                            0.376    
## absdiff.sqrtage+dat$shift.mar*(male==0)   <1e-04 ***
## nodefactor.deg.other.binary.1             <1e-04 ***
## offset(nodematch.male)                    <1e-04 ***
## offset(nodefactor.olderpartnerM.1)        <1e-04 ***
## offset(concurrent)                        <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
##  The following terms are fixed by offset and are not estimated:
##   offset(netsize.adj) offset(nodematch.male) offset(nodefactor.olderpartnerM.1) offset(concurrent)

4.2 MCMC Diagnostics / GOF

mcmc.diagnostics(fit.marriage)

gof <- gof(fit.marriage, GOF = "model")
saveRDS(gof, here("fits", "gof.m.rds"))
gof <- readRDS(here("fits", "gof.m.rds"))
plot(gof)

4.3 Static Diagnostics

static.mar <- netdx(marriage.netest, nsims = 500, dynamic = FALSE)
saveRDS(static.mar, here("fits", "static.m.rds"))
static.mar <- readRDS(here("fits", "static.m.rds"))
static.mar
## EpiModel Network Diagnostics
## =======================
## Diagnostic Method: Static
## Simulations: 500
## 
## Formation Diagnostics
## ----------------------- 
##                                              Target    Sim Mean Pct Diff
## edges                                      3478.196    3474.476   -0.107
## nodecov.ageF                             115489.827  115372.034   -0.102
## nodecov.ageF^2                          3960149.254 3956455.334   -0.093
## nodecov.ageF<30                            4512.671    4506.556   -0.136
## nodecov.ageF>30                            2282.498    2280.464   -0.089
## nodecov.ageM                             120933.300  120815.776   -0.097
## nodecov.ageM^2                          4326955.505 4322977.556   -0.092
## nodecov.ageM<30                            4216.771    4211.118   -0.134
## nodecov.ageM>30                            2562.346    2561.888   -0.018
## absdiff.sqrtage+dat$shift.mar*(male==0)     819.197     820.059    0.105
## nodefactor.deg.other.binary.1                24.909      24.984    0.303
## offset(nodematch.male)                           NA       0.000       NA
## offset(nodefactor.olderpartnerM.1)               NA       0.000       NA
## offset(concurrent)                               NA       0.000       NA
##                                            Sim SD
## edges                                      29.796
## nodecov.ageF                              969.749
## nodecov.ageF^2                          34934.018
## nodecov.ageF<30                            46.456
## nodecov.ageF>30                            23.726
## nodecov.ageM                             1001.210
## nodecov.ageM^2                          36557.584
## nodecov.ageM<30                            45.093
## nodecov.ageM>30                            25.112
## absdiff.sqrtage+dat$shift.mar*(male==0)    15.012
## nodefactor.deg.other.binary.1               4.965
## offset(nodematch.male)                      0.000
## offset(nodefactor.olderpartnerM.1)          0.000
## offset(concurrent)                          0.000
plot(static.mar)

4.4 Dynamic Diagnostics

dynamic.marriage <- netdx(marriage.netest, nsims = 5, nsteps = 1000, ncores = np, 
                       set.control.ergm = control.simulate.ergm(MCMC.burnin=1e5, 
                                                                MCMC.interval=1e5),
                       set.control.stergm = control.simulate.network(MCMC.burnin.min=1e5,
                                                                     MCMC.burnin.max=1e5))
saveRDS(dynamic.marriage, here("fits", "dynamic.m.rds"))
dynamic.marriage <- readRDS(here("fits", "dynamic.m.rds"))
plot(dynamic.marriage)
plot(dynamic.marriage)

plot(dynamic.marriage, type = "duration")

plot(dynamic.marriage, type = "dissolution")

4.5 Checking Mean Deg Dist by Age/Sex and Age/Sex Mixing are close to data

n <- dynamic.marriage$nw

plot(summary(n ~ nodefactor(~ageM)), type="l", main="Males")
lines(m, pch=19, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)

plot(summary(n ~ nodefactor(~ageF)), type="l", main="Females")
lines(f, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)

# age/sex mixing
aF <- n %v% "ageF"
aM <- n %v% "ageM" 

el <- as.edgelist(n)
fem1 <- aF[el[,1]]
fem2 <- aF[el[,2]]
male1 <- aM[el[,1]]
male2 <- aM[el[,2]]

allFs <- NULL
allMs <- NULL 

for (i in 1:length(fem1)){
  if (fem1[i] != 0) {allFs[i] <- fem1[i]} 
  else {allFs[i] <- fem2[i]}
  
  if (male1[i] != 0) {allMs[i] <- male1[i]} 
  else {allMs[i] <- male2[i]}
  }

agesel <- as.data.frame(cbind(allFs, allMs))
sqrtages <- sqrt(agesel)
sqrtages$diff <- sqrtages[,1] - sqrtages[,2]
mean(sqrtages$diff)
## [1] -0.1319895

5 Cohab Network

5.1 Ergm.Ego Estimation & Setup

x <- control.ergm.ego(ppopsize = 1, ppop.wt = "sample",
                      ergm.control = control.ergm(MCMLE.maxit = 200,
                                                  MCMC.interval = MCMCinterval,
                                                  MCMC.burnin = MCMCburnin))

x$ppopsize <- mpop

fit.cohab <- ergm.ego(egodat_cohab ~ edges + 
                         nodecov(~ageF) + 
                         nodecov(~ageF^2) +
                         nodecov(~ageM) + 
                         nodecov(~ageM^2) +
                         nodefactor(~ageF, levels=5:9) + 
                         nodefactor(~ageM, levels=7:10) + 
                         absdiff(~sqrtage + dat$shift.cohab*(male==0)) + 
                         nodefactor("deg.other.binary") +
                         offset(nodematch("male", diff = FALSE)) +
                         offset(nodefactor("olderpartnerC")) +
                         #offset(nodefactor("debuted", levels=1)) + 
                         offset("concurrent"), 
                         offset.coef = c(-Inf, -Inf, -Inf), #offset.coef = c(-Inf, -Inf, -Inf, -Inf),
                         control = x)

saveRDS(fit.cohab, here("fits", "fit.cohab.rds"))
fit.cohab <- readRDS(here("fits", "fit.cohab.rds"))
cohab.netest <- ego.netest(fit.cohab, DissolutionCohab)
summary(fit.cohab)
## Call:
## ergm(formula = ergm.formula, offset.coef = ergm.offset.coef, 
##     target.stats = m, eval.loglik = FALSE, control = control$ergm.control)
## 
## Iterations:  2 out of 200 
## 
## Monte Carlo MLE Results:
##                                             Estimate Std. Error MCMC % z value
## offset(netsize.adj)                       -9.9034876  0.0000000      0    -Inf
## edges                                     -8.7874426  0.9013938      0  -9.749
## nodecov.ageF                               0.4025384  0.0667472      0   6.031
## nodecov.ageF^2                            -0.0077624  0.0012306      0  -6.308
## nodecov.ageM                               0.3009793  0.0578495      0   5.203
## nodecov.ageM^2                            -0.0046020  0.0009947      0  -4.627
## nodefactor.ageF.18                         0.1432191  0.3019090      0   0.474
## nodefactor.ageF.19                         0.6184190  0.2804912      0   2.205
## nodefactor.ageF.20                         1.5189275  0.3457406      0   4.393
## nodefactor.ageF.21                         0.8620661  0.2534431      0   3.401
## nodefactor.ageF.22                         0.6722736  0.2334703      0   2.879
## nodefactor.ageM.20                         0.5421379  0.3020196      0   1.795
## nodefactor.ageM.21                         0.3026223  0.2997873      0   1.009
## nodefactor.ageM.22                         0.7057102  0.2506985      0   2.815
## nodefactor.ageM.23                         0.8937318  0.2602071      0   3.435
## absdiff.sqrtage+dat$shift.cohab*(male==0) -2.7147838  0.1335418      0 -20.329
## nodefactor.deg.other.binary.1             -2.4469524  0.3307137      0  -7.399
## offset(nodematch.male)                          -Inf  0.0000000      0    -Inf
## offset(nodefactor.olderpartnerC.1)              -Inf  0.0000000      0    -Inf
## offset(concurrent)                              -Inf  0.0000000      0    -Inf
##                                           Pr(>|z|)    
## offset(netsize.adj)                        < 1e-04 ***
## edges                                      < 1e-04 ***
## nodecov.ageF                               < 1e-04 ***
## nodecov.ageF^2                             < 1e-04 ***
## nodecov.ageM                               < 1e-04 ***
## nodecov.ageM^2                             < 1e-04 ***
## nodefactor.ageF.18                        0.635230    
## nodefactor.ageF.19                        0.027470 *  
## nodefactor.ageF.20                         < 1e-04 ***
## nodefactor.ageF.21                        0.000670 ***
## nodefactor.ageF.22                        0.003983 ** 
## nodefactor.ageM.20                        0.072647 .  
## nodefactor.ageM.21                        0.312756    
## nodefactor.ageM.22                        0.004878 ** 
## nodefactor.ageM.23                        0.000593 ***
## absdiff.sqrtage+dat$shift.cohab*(male==0)  < 1e-04 ***
## nodefactor.deg.other.binary.1              < 1e-04 ***
## offset(nodematch.male)                     < 1e-04 ***
## offset(nodefactor.olderpartnerC.1)         < 1e-04 ***
## offset(concurrent)                         < 1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
##  The following terms are fixed by offset and are not estimated:
##   offset(netsize.adj) offset(nodematch.male) offset(nodefactor.olderpartnerC.1) offset(concurrent)

5.2 MCMC Diagnostics / GOF

mcmc.diagnostics(fit.cohab)

gof <- gof(fit.cohab, GOF = "model")
saveRDS(gof, here("fits", "gof.c.rds"))
gof <- readRDS(here("fits", "gof.c.rds"))
plot(gof)

5.3 Static Diagnostics

static.cohab <- netdx(cohab.netest, nsims = 500, dynamic = FALSE)
saveRDS(static.cohab, here("fits", "static.c.rds"))
static.cohab <- readRDS(here("fits", "static.c.rds"))
static.cohab
## EpiModel Network Diagnostics
## =======================
## Diagnostic Method: Static
## Simulations: 500
## 
## Formation Diagnostics
## ----------------------- 
##                                                Target    Sim Mean Pct Diff
## edges                                        1334.405    1333.460   -0.071
## nodecov.ageF                                36428.502   36382.962   -0.125
## nodecov.ageF^2                            1048047.957 1046062.370   -0.189
## nodecov.ageM                                39792.620   39756.850   -0.090
## nodecov.ageM^2                            1247690.001 1246314.762   -0.110
## nodefactor.ageF.18                             26.033      25.994   -0.149
## nodefactor.ageF.19                             58.637      58.818    0.309
## nodefactor.ageF.20                             81.757      81.176   -0.711
## nodefactor.ageF.21                             88.739      89.436    0.786
## nodefactor.ageF.22                             83.736      84.506    0.919
## nodefactor.ageM.20                             40.588      40.286   -0.745
## nodefactor.ageM.21                             36.329      36.610    0.774
## nodefactor.ageM.22                             89.100      89.390    0.325
## nodefactor.ageM.23                             76.466      76.920    0.594
## absdiff.sqrtage+dat$shift.cohab*(male==0)     418.638     417.453   -0.283
## nodefactor.deg.other.binary.1                  60.629      60.448   -0.298
## offset(nodematch.male)                             NA       0.000       NA
## offset(nodefactor.olderpartnerC.1)                 NA       0.000       NA
## offset(concurrent)                                 NA       0.000       NA
##                                              Sim SD
## edges                                        30.870
## nodecov.ageF                                872.239
## nodecov.ageF^2                            27193.786
## nodecov.ageM                                948.242
## nodecov.ageM^2                            31918.371
## nodefactor.ageF.18                            4.686
## nodefactor.ageF.19                            6.923
## nodefactor.ageF.20                            6.859
## nodefactor.ageF.21                            8.038
## nodefactor.ageF.22                            7.621
## nodefactor.ageM.20                            5.809
## nodefactor.ageM.21                            5.477
## nodefactor.ageM.22                            8.211
## nodefactor.ageM.23                            7.485
## absdiff.sqrtage+dat$shift.cohab*(male==0)    14.249
## nodefactor.deg.other.binary.1                 7.512
## offset(nodematch.male)                        0.000
## offset(nodefactor.olderpartnerC.1)            0.000
## offset(concurrent)                            0.000
plot(static.cohab)

5.4 Dynamic Diagnostics

dynamic.cohab <- netdx(cohab.netest, nsims = 5, nsteps = 1000, ncores = np, 
                       set.control.ergm = control.simulate.ergm(MCMC.burnin=1e5, 
                                                                MCMC.interval=1e5),
                       set.control.stergm = control.simulate.network(MCMC.burnin.min=1e5,
                                                                     MCMC.burnin.max=1e5))
saveRDS(dynamic.cohab, here("fits", "dynamic.c.rds"))
dynamic.cohab <- readRDS(here("fits", "dynamic.c.rds"))
plot(dynamic.cohab)
plot(dynamic.cohab)

plot(dynamic.cohab, type = "duration")

plot(dynamic.cohab, type = "dissolution")

5.5 Checking Mean Deg Dist by Age/Sex and Age/Sex Mixing are close to data

n <- dynamic.cohab$nw

plot(summary(n ~ nodefactor(~ageM)), type="l", main="Males")
lines(mc, pch=19, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)

plot(summary(n ~ nodefactor(~ageF)), type="l", main="Females")
lines(fc, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)

# age/sex mixing
aF <- n %v% "ageF"
aM <- n %v% "ageM" 

el <- as.edgelist(n)
fem1 <- aF[el[,1]]
fem2 <- aF[el[,2]]
male1 <- aM[el[,1]]
male2 <- aM[el[,2]]

allFs <- NULL
allMs <- NULL 

for (i in 1:length(fem1)){
  if (fem1[i] != 0) {allFs[i] <- fem1[i]} 
  else {allFs[i] <- fem2[i]}
  
  if (male1[i] != 0) {allMs[i] <- male1[i]} 
  else {allMs[i] <- male2[i]}
  }

agesel <- as.data.frame(cbind(allFs, allMs))
sqrtages <- sqrt(agesel)
sqrtages$diff <- sqrtages[,1] - sqrtages[,2]
mean(sqrtages$diff)
## [1] -0.2344093

6 Other Network

6.1 Fit

fit.other <- ergm.ego(egodat_other ~ edges + 
                         nodecov(~ageF) + 
                         nodecov(~ageF<19) + 
                         nodecov(~ageF^2) +
                         nodecov(~ageM<19) + 
                         nodecov(~ageM) + 
                         nodecov(~ageM^2) +
                        absdiff(~sqrtage + dat$shift.other*(male==0)) + 
                        nodefactor("deg.mar") +
                        nodefactor("deg.cohab") +
                        concurrent(by=~male) +
                        offset(nodematch("male", diff = FALSE)) +
                        offset(nodefactor("olderpartnerO")),
                        offset.coef = c(-Inf, -Inf), 
                        control = x)


#offset(nodefactor("debuted", levels=1)),
#offset.coef = c(-Inf, -Inf),

saveRDS(fit.other, here("fits", "fit.other.rds"))
fit.other <- readRDS(here("fits", "fit.other.rds"))
other.netest <- ego.netest(fit.other, DissolutionOther)
summary(fit.other)
## Call:
## ergm(formula = ergm.formula, offset.coef = ergm.offset.coef, 
##     target.stats = m, eval.loglik = FALSE, control = control$ergm.control)
## 
## Iterations:  3 out of 200 
## 
## Monte Carlo MLE Results:
##                                            Estimate Std. Error MCMC % z value
## offset(netsize.adj)                       -9.903488   0.000000      0    -Inf
## edges                                     -6.065217   1.606071      0  -3.776
## nodecov.ageF                               0.448504   0.083211      0   5.390
## nodecov.ageF<19                           -0.337883   0.228428      0  -1.479
## nodecov.ageF^2                            -0.008739   0.001446      0  -6.045
## nodecov.ageM<19                           -0.050555   0.200336      0  -0.252
## nodecov.ageM                               0.270125   0.074329      0   3.634
## nodecov.ageM^2                            -0.003829   0.001247      0  -3.071
## absdiff.sqrtage+dat$shift.other*(male==0) -3.103887   0.128697      0 -24.118
## nodefactor.deg.mar.1                      -4.387958   0.335242      0 -13.089
## nodefactor.deg.cohab.1                    -3.673664   0.330056      0 -11.130
## concurrent.male0                          -3.374662   0.272150      0 -12.400
## concurrent.male1                          -1.471438   0.299298      0  -4.916
## offset(nodematch.male)                         -Inf   0.000000      0    -Inf
## offset(nodefactor.olderpartnerO.1)             -Inf   0.000000      0    -Inf
##                                           Pr(>|z|)    
## offset(netsize.adj)                        < 1e-04 ***
## edges                                     0.000159 ***
## nodecov.ageF                               < 1e-04 ***
## nodecov.ageF<19                           0.139096    
## nodecov.ageF^2                             < 1e-04 ***
## nodecov.ageM<19                           0.800770    
## nodecov.ageM                              0.000279 ***
## nodecov.ageM^2                            0.002132 ** 
## absdiff.sqrtage+dat$shift.other*(male==0)  < 1e-04 ***
## nodefactor.deg.mar.1                       < 1e-04 ***
## nodefactor.deg.cohab.1                     < 1e-04 ***
## concurrent.male0                           < 1e-04 ***
## concurrent.male1                           < 1e-04 ***
## offset(nodematch.male)                     < 1e-04 ***
## offset(nodefactor.olderpartnerO.1)         < 1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
##  The following terms are fixed by offset and are not estimated:
##   offset(netsize.adj) offset(nodematch.male) offset(nodefactor.olderpartnerO.1)

6.2 MCMC Diagnostics / GOF

mcmc.diagnostics(fit.other)

gof.o <- gof(fit.other, GOF = "model")
saveRDS(gof.o, here("fits", "gof.o.rds"))
gof.o <- readRDS(here("fits", "gof.o.rds"))
plot(gof.o)

6.3 Static Diagnostic

static.other <- netdx(other.netest, nsims = 500, dynamic = FALSE)
## 
## Network Diagnostics
## -----------------------
## - Simulating 500 networks
## - Calculating formation statistics
saveRDS(static.other, here("fits", "static.other.rds"))
static.other <- readRDS(here("fits", "static.other.rds"))
static.other
plot(static.other)

6.4 Dynamic Diagnostic

dynamic.other <- netdx(other.netest, nsims = 5, nsteps = 1000, ncores = np,
                       set.control.ergm = control.simulate.ergm(MCMC.burnin=1e5, 
                                                                MCMC.interval=1e5),
                       set.control.stergm = control.simulate.network(MCMC.burnin.min=1e5,
                                                                     MCMC.burnin.max=1e5))
saveRDS(dynamic.other, here("fits", "dynamic.other.rds"))
dynamic.other <- readRDS(here("fits", "dynamic.other.rds"))
dynamic.other
## EpiModel Network Diagnostics
## =======================
## Diagnostic Method: Dynamic
## Simulations: 5
## Time Steps per Sim: 1000
## 
## Formation Diagnostics
## ----------------------- 
##                                                Target    Sim Mean Pct Diff
## edges                                        1951.961    1949.035   -0.150
## nodecov.ageF                                48913.199   48823.057   -0.184
## nodecov.ageF<19                              2287.271    2285.052   -0.097
## nodecov.ageF^2                            1321309.489 1318712.896   -0.197
## nodecov.ageM<19                              2181.387    2180.544   -0.039
## nodecov.ageM                                52313.794   52230.803   -0.159
## nodecov.ageM^2                            1511564.771 1509422.882   -0.142
## absdiff.sqrtage+dat$shift.other*(male==0)     553.355     553.436    0.015
## nodefactor.deg.mar.1                           60.275      60.282    0.012
## nodefactor.deg.cohab.1                         63.409      66.256    4.489
## concurrent.male0                               41.854      40.940   -2.184
## concurrent.male1                              146.741     147.834    0.745
## offset(nodematch.male)                             NA       0.000       NA
## offset(nodefactor.olderpartnerO.1)                 NA       0.000       NA
##                                              Sim SD
## edges                                        33.197
## nodecov.ageF                                879.840
## nodecov.ageF<19                              41.388
## nodecov.ageF^2                            26937.054
## nodecov.ageM<19                              40.336
## nodecov.ageM                                908.185
## nodecov.ageM^2                            29079.782
## absdiff.sqrtage+dat$shift.other*(male==0)    13.786
## nodefactor.deg.mar.1                          7.517
## nodefactor.deg.cohab.1                        7.645
## concurrent.male0                              6.279
## concurrent.male1                             11.095
## offset(nodematch.male)                        0.000
## offset(nodefactor.olderpartnerO.1)            0.000
## 
## Dissolution Diagnostics
## ----------------------- 
##                Target Sim Mean Pct Diff Sim SD
## Edge Duration  55.000   51.781   -5.853 51.175
## Pct Edges Diss  0.018    0.018    0.167  0.003
plot(dynamic.other)

plot(dynamic.other, type = "duration")

plot(dynamic.other, type = "dissolution")

6.5 Checking Mean Deg Dist by Age/Sex and Age/Sex Mixing are close to data

n <- dynamic.other$nw

# mean deg dist 
plot(summary(n ~ nodefactor(~ageM)), type="l", main="Males")
lines(mo, pch=19, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)

plot(summary(n ~ nodefactor(~ageF)), type="l", main="Females")
lines(fo, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)

# age/sex mixing
aF <- n %v% "ageF"
aM <- n %v% "ageM" 

el <- as.edgelist(n)
fem1 <- aF[el[,1]]
fem2 <- aF[el[,2]]
male1 <- aM[el[,1]]
male2 <- aM[el[,2]]

allFs <- NULL
allMs <- NULL 

for (i in 1:length(fem1)){
  if (fem1[i] != 0) {allFs[i] <- fem1[i]} 
  else {allFs[i] <- fem2[i]}
  
  if (male1[i] != 0) {allMs[i] <- male1[i]} 
  else {allMs[i] <- male2[i]}
  }

agesel <- as.data.frame(cbind(allFs, allMs))
sqrtages <- sqrt(agesel)
sqrtages$diff <- sqrtages[,1] - sqrtages[,2]
mean(sqrtages$diff)
## [1] -0.1753593
#netests <- list(marcoh.netest, other.netest)
#saveRDS(netests, here("Setup", "FinalNetworks", "objects", "full", "nodefactor", "netests_nodefactor.rds"))